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Abstract A birth-death process is a continuous-time Markov chain that counts the 
number of particles in a system over time. In the general process with n current par- 
ticles, a new particle is born with instantaneous rate \ n and a particle dies with 
instantaneous rate fin- Currently no robust and efficient method exists to evaluate 
the finite-time transition probabilities in a general birth-death process with arbitrary 
birth and death rates. In this paper, we first revisit the theory of continued fractions 
to obtain expressions for the Laplace transforms of these transition probabilities and 
make explicit an important derivation connecting transition probabilities and contin- 
ued fractions. We then develop an efficient algorithm for computing these probabilities 
that analyzes the error associated with approximations in the method. We demonstrate 
that this error-controlled method agrees with known solutions and outperforms previ- 
ous approaches to computing these probabilities. Finally, we apply our novel method 
to several important problems in ecology, evolution, and genetics. 
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1 Introduction 

Birth-death processes (BDPs) have a rich history in probabilistic modeling, includ- 
ing applications in ecology, genetics, and evolution ( jThorne et al 1991 Krone and 
Neuhauser|1997 Novozhilov et al|2006 |. Traditionally, BDPs have been used to model 
the number of organisms or particles in a system, each of which reproduce and die in 
continuous time. A general BDP is a continuous-time Markov chain on the non-negative 
integers in which instantaneous transitions from state n > to either n + 1 or n — 1 are 
possible. These transitions are called "births" and "deaths". Starting at state n, jumps 
to n + 1 occur with instantaneous rate X n and jumps to n — 1 with instantaneous 
rate fi n - The simplest BDP has linear rates An = nX and fi n = n\i with no state- 



independent terms (Kendall 1948 Feller 19711. This model is the most widely- used 



BDP since there exist closed-form expressions for its transition probabilities (Bailey 



1964] |Novozhilov et al|2006l >. Many applications of BDPs require convenient methods 
for computing the probability P m ,n(t) that the system moves from state m to state n 
in finite time t > 0. These probabilities exhibit their usefulness in many modeling ap- 
plications since the probabilities do not depend on the possibly unobserved path taken 
by the process from m to n and hence make possible analyses of discretely sampled or 
partially observed processes. Despite the relative simplicity of specifying the rates of a 
general BDP, it can be remarkably difficult to find closed-form solutions for the tran- 
sition probabilities even for simple models (Rcnshaw 1993; Mcdcrcr 2003, Novozhilov 
et al|2006 l. 

In a pioneering series of papers, Karlin a nd McGregor| develop a formal theory 
of general BDPs that expresses their transition probabilities in terms of a sequence of 
orthogonal polynomials and a spectral measure (Karlin and McGregor|1957a|b|[l958bl ). 
While the work of Karlin and McGregor yields valuable theoretical insights regarding 
the existence of unique solutions and properties of recurrence and transience for a given 
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process, there remains no clear recipe for determining the orthogonal polynomials and 
measure corresponding to an arbitrary set of birth and death rates. Additionally, even 
when the polynomials and measure are known, the transition probabilities may not 
have an analytic representation or a convenient computational form. 

Possibly due to the difficulty of finding computationally useful formulas for transi- 
tion probabilities in general BDPs, many applied researchers resort to easier analyses 
using moments, first passage times, equilibrium probabilities, and other tractable quan- 
tities of interest. Referring to the system of Kolmogorov forward differential equations 



for transition probabilities that we give below, Novozhilov et al ( 2006 page 73) write 



"The problem with exact solutions of system (1) is that, in many cases, the 
expressions for the state probabilities, although explicit, are intractable for 
analysis and include special polynomials. In such cases, it may be sensible 
to solve more modest problems concerning the birth-and-death process under 
consideration, without the knowledge of the time-dependent behavior of state 
probabilities p n (t)." 

Indeed, closed-form analytic expressions for transition probabilities of general BDPs 
are only known for a few types of processes. Some examples include constant birth and 
death rates (Bailey 19641, zero birth or death rates (pure-death and pure-birth) (Yule 



1925 Taylor and Karlin|1998 l, and certain linear rates ( Karlin and McGregor|1958a l 



As a seemingly straightforward example, in the BDP with linear birth and death rates 



A n = n\ + v and fi n = rifJ, + 7 including state-independent terms, Ismail et al ( 1988) 
offer the orthogonal polynomials and associated measure, but still no closed form is 
available for the transition probabilities. 

Despite the difficulty in obtaining analytic expressions, several authors have made 
progress in approximate numerical methods for solution of transition probabilities 



in general BDPs. Murphy and O'Donohoe ( 1975 1 develop an appealing numerical 



method for the transition probabilities based on a continued fraction representation 
of Laplace-transformed transition probabilities. They invert these transformed prob- 
abilities by first truncating the continued fraction. Several other authors give similar 
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expressions derived from truncation of the state space ( |Grassmann|1977b|a Rosenlund 



p78] |Sharma and Dass|p88] |Mohanty et~al|p93| |. However, |Klar et al| p010| ) find 
that methods based on continued fraction truncation and then subsequent analytical 
transformation can suffer from instability. As an alternative, [Parthasarathy lind Sud-] 



hesh (2006a I express the infinite continued fraction representation given by Murphy 
|and O'Donohoe] as a power series. Unfortunately, the small radius of convergence of 
this series makes it less useful for numerical computation. 

We also note that for general BDPs that take values on a finite state space (usually 
n 6 {0, f , . . . , N}), it is possible to write a finite-dimensional stochastic transition rate 
matrix and solve for the matrix of transition probabilities. If the rate matrix is diagonal- 
izable, computation of transition probabilities in this manner can be computationally 
straightforward. To illustrate, let Q be a finite-dimensional stochastic rate matrix with 
Q = UAU" 1 where U is an orthogonal matrix and A is diagonal. The matrix of tran- 
sition probabilities P satisfies the matrix differential equation P = PQ with initial 
condition P(0) = I. The solution is P(t) = exp[Qt] = U diag(e Zl ', e Z2 \ e Zjv *) U~ X , 
where Z\, . . . , zjv are the eigenvalues of Q. However, it is possible to specify reason- 
able rate parameters in a general BDP that satisfy requirements for the existence of 
a unique solution, but do not result in a diagonalizable rate matrix. Also, if the state 
space over which the BDP takes values is large, numerical eigendecomposition of Q 
may be computationally expensive and could introduce serious roundoff errors. 

To our knowledge, no robust computational method currently exists for finding the 
finite-time transition probabilities of general BDPs with arbitrary rates. Such a tech- 
nique would allow rapid development of rich and sophisticated ecological, genetic, and 
evolutionary models. Additionally, in statistical applications, transition probabilities 
can serve as observed data likelihoods, and are thus often useful in estimating tran- 
sition rate parameters from partially observed BDPs. We believe more sophisticated 
BDPs can be very useful for applied researchers. In spite of the numerical difficulties 
presented by approximant methods, we are surprised that continued fraction meth- 



ods like that of Murphy and O'Donohoe ( 1975 1 are not more widely explored. This 
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may be due to omission of important details in their derivation of continued fraction 
expressions for the Laplace transform of the transition probabilities. 

In this paper, we build on continued fraction expressions for the Laplace transforms 
of the transition probabilities of a general BDP using techniques similar to those in- 
troduced by Murphy and O'Donohoe, and we fill in the missing details in the proof of 
this representation. We then apply the Laplace inversion formulae of |Abate and Whitt| 
(1992a b) to obtain an efficient and robust method for computation of transition prob- 
abilities in general BDPs. Our method relies on three observations: 1) it is possible 
to find exact expressions for Laplace transforms of the transition probabilities of a 
general BDP using continued fractions ( jMurphy and O'Donoh oc 1975); 2) evaluation 
of continued fractions is typically very fast, requires far fewer evaluations than equiv- 
alent power series, and there exist robust algorithms for evaluating them efficiently 



( Bankier and Leighton|l942"| |Wall|1 948 Bla nch|1964||Lorentzen and Waadeland|1992 



Craviotto et al||1993| |Abate and Whitt|[l999| |Cuyt et al||2008| ; and 3) recovery of 
probability distributions by Laplace inversion using a Riemann sum approximation is 



often more computationally stable than analytical methods of inversion (Abate and 



Whitt 1992a b 19951. Finally, we demonstrate the advantages of our error-controlled 
method through its application to several birth-death models in ecology, genetics, and 
evolution whose solution remains unavailable by other means. 



2 Transition probabilities 

2.1 Background 

A general birth-death process is a continuous-time Markov process X — {X(t),t > 0} 
counting the number of arbitrarily defined "particles" in existence at time t > 0, with 
X(0) — m > 0. To characterize the process, we define non-negative instantaneous 
birth rates A„ and death rates \i n for n > 0, with fiQ — and transition probabilities 
Pm,n(t) = P r (X(t) = n \ X(0) — m). While X n and /j, n are time- homogeneous con- 
stants, they may depend on n. We refer to the classical linear BDP in which \ n = nX 
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and fi n = rifx as the "simple birth-death process" (Kendall 1948 Feller| 1971 1. The 
general BDP transition probabilities satisfy the infinite system of ordinary differential 
equations 



dPm.oft) 

dt 

dP m ,n(t) 

dt 



fl-P 71,1 

(t) - Ao-P m ,o(t), and 
>>n-lPm,n-l(t) + Hn+lPm.n+lif) - (An + Hn)Pm,n{t) for n > 1, 



(1) 



with boundary conditions -Pm,m(0) — 1 and Pm,n (0) = for n / m (|Feller|p7T|. 



Karlin and McGregor ( 1957b I show that for arbitrary starting state m, transition 



probabilities can be represented in the form 

f-OO 

Pm,n(t) = 7r„ / e~ xt Q m (x)Q n {x)iP{dx), (2) 
JO 

where 7ro = 1 and 7r n = (Aq • ■ ■ \ n —l)/ (m ■ ' ■ A 4 ") f° r m > 1. Here, {Qn(x)} is a sequence 
of polynomials satisfying the three-term recurrence relation 

^oQl(x) = Ao + /io — a;, and 

(3) 

AnQn+l(l) = (An + £in — x)Q n (x) ~ fl n Qn-l{x), 

and ^) is the spectral measure of X with respect to which the polynomials {Qn(x)} are 
orthogonal. The system |l]) has a unique solution if and only if 

E («* + A^) = °°- (4) 



fe=0 



In what follows, we assume that the rate parameters {An} and {fin} satisfy Q. Closed- 
form solutions to |l]) are available for a surprisingly small number of choices of {An} 
and {fin}- We therefore need another approach to find useful formulae for computation 
of the transition probabilities. 
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2.2 Continued fraction representation of Laplace transform 

To find an expression that is useful for computing P m ,n(t) for an arbitrary general BDP, 
a fruitful approach is often to Laplace transform each equation of the system and 
form a recurrence relationship relating back to the Laplace transform of P m ,n(t)- We 



base our presentation on that of Murphy and O'Donohoe (19751. Denote the Laplace 
transform of P n ,m(t) as 



fm,n( S ) = £ [Pm.nit)] 00 = / e P m ,„(t) At. (5) 

JO 

Applying the Laplace transform to |TJ, with the starting state m = 0, we arrive at 



s/o,o(s) - fo,o(0) = /til/o.lOO ~ Wo.oOOi and 

sfo,n{s) - fb,n(0) = A„_i/ 0: „-l(s) + Mri + l/0,n+l(s) - (A„ + fln)/o,ti(s) 



(6) 



for n > 1. Rearranging and recalling that Po,o(0) = 1 and Po,n(0) = for n > 1, we 
simplify Q to 



/o,l00 = — [O 5 + Ao)/o,o(s) - l],and 



/t),n00 



1 



(s + A„_i + Atri-l)/o,n-l(s) - A„_ 2 /o,ri-2(s) 



for n > 2. 



Some rearranging of |7|) yields the forward system of recurrence relations 

1 



/o,o00 

/o,n00 



A„,_l 



and 



(7) 



(8) 



Then combining these expressions, we arrive at the generalized continued fraction 



/o,o00 = 



(9) 



s + Ao - 



Aqmi 



S + Al + /Ltl 



s + A 2 + ^2 — 



s 



This is an exact expression for the Laplace transform of the transition probability 
Po,o(t)- Let the partial numerators in Q be a\ — 1 and a n — — A n _2/i n — 1, an d the 
partial denominators 61 = s + \q and b n — s + \ n —i + fJ, n —l f° r n > 2. Then (JoJ 
becomes 

"1 



/o,o(s) 



61 + 



a 2 



62 



"3 
&3 + ' 



To express ( 10 1 in more typographically economical notation, we write 



(10) 



/o,o(s) = 



ai a-2 a-j 
h + 62+ 63+ 



(11) 



We denote the fcth convergent (approximant) of /o,o( s ) as 



'0,0 W~ bl+b2+ 



(12) 



There are deep connections between the orthogonal polynomial representation 
Laplace transforms and continued fractions of the form ||9| that are beyond the 
scope of this paper ( [Karlin and McGre gor 1957b Bordes and Roehner|19 83 Guillemin 
|and Pinchon|[l9 99;). Interestingly, [FTajolet and Guillemin| ( |2000[ ) demonstrate a close 
relationship between the Laplace transforms of transition probabilities and state paths 
of the underlying Markov chain. 

Before stating a theorem supporting this representation, we give two lemmas that 
will be useful in what follows. 



Lemma 1 Both the numerator and denominator of ( 12 1 satisfy the same re 



currence, due to Wallis 1695) 



A k = b k A k-l + a k A k-2, and 
Bk = bkBk-i + ctfc-Bfe-2, 



(13) 



with Aq = 0, Ai — a\, Bq = 1, and B\ = b\. 
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Lemma 2 By repeated application of Lemrna^ we arrive at the determinant formula 



AkBk-i — Af.-iBf. = {b k A k _ 1 + a k A k _ 2 )Bk-i — 4fc-i(&fc-Bfc-i + a k B k _ 2 ) 
= — a k{A k -\B k _2 — -Afc-a-Bfc-l) 
= (-l) k - 1 f[a i . 



(14) 



Now we state and prove a theorem giving expressions for the Laplace transform of 



Pm,n{t). Although Murphy and O'Donohoe (19751 first report this result, they do not 
provide a detailed derivation in their paper. 

Theorem 1 The Laplace transform of the transition probability P m ,n(t) is given by 

B n (s) B m (s)a m+ 2 a m+ 3 



fm,n{s) — * 



n 

\j=n+l 



B m +i(s)+ b m+ 2+ &m+3+ 



B m (s) B n (s)a n+2 a n+ 3 



I n-l ' 

n X i J B n+1 (s)+ b n+2 + b n+3 + 



for n <m, 



for m < n, 



(15) 



where a n , b n , and B n are as defined above. 

Proof To simplify notation, we sometimes omit the dependence of f k , A k , and B k on 
the Laplace variable s. Suppose the process starts at X(0) — m. We can re-write the 
Laplace-transformed equations (JsJ with P m ,m(0) = 1 and P m ,n(0) = for all n 7^ m 
as 

sfm,o( s ) ~ S m0 = Ml/m,l( s ) - Wro,o(s), (16a) 
sfm,n(s) - 8mn = Ki-lfm,n-l (s) + A*n+l/m,n+l( s ) ~~ C^n + ^n)fm.n{s), (16b) 



where <5 mn = 1 if m = n and zero otherwise. We first derive the expression for n < m. 



If m = 0, /o,o( s ) is given by (111, so we assume in what follows when n < m, that 



m > 1. Rearranging ( 16a I, we see that since Bq — 1 and s + Aq = 61 = Si, 



r -BO r 

B\ 



(17) 
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Now, to show the general case by induction, assume that for n < m, 



f _ B n -i 

Jm.n — l — D A*n/m,7v 
tin 



(18) 



Substituting (181 into (16b I when n < m, we have 



tin 



+ a n + l~ g ) fm.n — Mn+l/m.n 



f — 71 f 
Jm.n — ^5 fn+l/m,Ti+l 

tin+1 



1+1 



and so ( 18 1 is true for any n < m. Letting n = m, we have by ( 18 1 and ( 16b l 



(19) 

(20) 
(21) 



bm+lfm.m — 1 + A m _ 1 ( — — j-imfm.m ) + Mrn+l/m,m+l • (22) 



Recalling that s + A m + /u m = b m +i and using LemmaJTJ 



r _ i -Bm+i f 

Pi7i+Um,ro+1 — 1 5 Jm,m- 
Dm 



(23) 



Rearranging the previous equation, we find that 



fm.m — 



Likewise, we can write (|16b| as a continued fraction recurrence: 



(24) 



fm.n -Vi — 1 

/m,ra-l S + /Xn + Art + Mrt+l J^"* 1 



(25) 



Then plugging (25 1 into (24 1 and iterating, we obtain the continued fraction for f m ,r. 



1 



1m+2 lm+3 



/ m , m — 



&m+2 + &771+3 + 



-Bm -Bmdm+2 Om+3 
-B m +1+ b m +2+ &m+3 + 



(26) 
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This is an exact formula for the Laplace transform of P m ,m(t), and proves the case 
m = n. For n < m, we iterate (|18[) to get 



fn 



t^n+1 /m,n+l 

Mn+lMri+2/m,n+2 

A t n+lMri+2 ' ■ ' l^mfn 



B n +1 

B n B n+ \ 
B n +i B n+ 2 
B n B n+ \ B m _i 



B n +i B n+ 2 B 
' " \ B, 



n w 



a. 



fn 



(27) 



Substituting (261 for f m ,m completes the proof for n<m. 



To find the formula for f m ,n when n > m, we adopt a similar approach. From (24 1 
we arrive at 

B m +lfm,m — B m — BmCm+l/m,m+l (28) 



We proceed inductively. Assume that for n > m, 



B n +lfm,n — I Y\ B m J r ^n+\Brifm,n+l- 



\3=m 



(29) 



From (16b I, we have 



b n +2fm,n+l — X n fm,n + Hn+2fm,n+2- 



(30) 



Solving for f m ,n in ( 29 1 and plugging this into the above equation, we have 



n-1 



Bn 
B n +1 



b n +2fm,n+l — An I J^J Xj | - — + XnHn+1 p ~ fm,n+l + A t ra+2/m,7t+2 ■ (31) 



Bn 
B n +1 ' 



Recalling that — Aji/i n -|_i = a n +2> 



{b n +2B n +l + a n +2B m ) fm,n+l = I Aj J £> n + fl n +2B n +lf; 



m,n+2 j 



(32) 
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and by Lemma [T] 

B n +2fm,n+l = I Aj B m + jtt m _|_2-B n +l fm,m+2 ■ (33) 



This establishes the recurrence ([29]). Then for any n > m, we can rearrange (|29j) to 
obtain 

fm,n = f if ~ /mn+1 - 04) 



This completes the proof. 



2.3 Obtaining transition probabilities 



Murphy and O'Donohoe ( 1975| find transition probabilities by truncating ( 15 1 at a pre- 
specified depth, forming a partial fractions sum, and inverse transforming. Parthasarathy 



and Sudhcsh (2006a) give a series solution for transition probabilities based on an equiv- 



alence between continued fractions like ( |15[ l and power series. However, both of these 
approaches suffer from serious drawbacks, as we explore in detail in the Appendix. 
We instead seek an efficient and robust numerical method for evaluating and in- 



verting (15). We first note that continued fractions typically converge rapidly, and in 



our experience, evaluation of ( 15 1 is very fast and stable using the Lentz algorithm and 



its subsequent improvements ( Lentzp976 Thompson and Barnctt 1986; Press|[2007 ) 
We therefore invert (1151) numerically by a summation formula. 



To do this, we treat the continued fraction representation ( 15 1 of the Laplace trans 



form of P m ,n(t) as an unknown but computable function of the complex Laplace vari- 



able s. We base our presentation on that of Abate and Whitt ( 1992a I. If e is a positive 



real number such that all singularities of f m ,n(s) lie to the left of e in the complex 
plane, the inverse Laplace transform of f m ,n(s) is given by the Bromwich integral 



l (t)=£~ 1 (/ m ,„(s)) = 



2« /- 



e+ioo 



e s fm,n(s) ds. 



(35) 
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Letting s = e + iu, 
1 f°° 

Pm,n(t) = — / e^ +m ^*/m,n( e + dtl 
27r ./-co 

— / [ cos(ut) + i sin(ui)l / m n (e + iu) du 
2n /-oo 

/ Re(/ m ,„(e + iu)) cos(ut) — Im(/ m , rl (e + iu)) sin(ut) du 

J — oo J 

+ i / [lm(/ m , n (£ + iu)) cos(ut) + Re(/ TOj „(e + iu)) sin(ut)J du 



e 

2^ 



(36) 



but Pm,n(t) is real- valued, so the imaginary part of the last equality in (36 1 is zero. 
Then 

P m ,n(t) = — y ^Re(/ m , n (e + iu)) cos(ui) - Im(/ m ,„(e + iu)) sin(ut)J du. (37) 
But since P m ,nif) = for t < 0, we also have that 

/ Re(/ m ,n(e + iu)) cos(ui) + Im(/ m ,„(e + iu)) sin(ui) du = 0. (38) 

J — oo 



Then applying (381 to (371, we obtain 



Pm,n(t) — 



Finally, we note that since 



Re 



— / Re(/ m , n (e + iu)) cos(ut) du. 

^ J — oo 



(39) 



POO 

}(/(e - iu)) = / e~ et cos(ut)P m ,„(t) d* = Ee(/(e + iu)) , 
•/o 

it must be the case that Re(/ m , n (e + iu)) is even in u for every e. Therefore, 

2e ct Z -00 

P m ,n(t) = / Re(/ m ,n(e + iu)) cos(ut) du. 

71 Jo 



(40) 



(41) 
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Following Abate and Whitt (1992a]), we approximate the integral above by a dis- 
crete Riemann sum via the trapezoidal rule with step size h: 



he et 2he et °° 
Pm,n(*) » Re (/m.n(e)) + V Re (fm,n(e + ikh)) cos(kht) 

7T 7T ^ — ' 

fe=l 



=^4/2 



2f 



-Re /, 



t 



£(-l) fc Re(/, 



fc=i 



A + 2kni 
2t 



(42) 



where the second line is obtained by setting h — rr/(2t) and e = A/(2t); this change of 
variables eliminates the cosine term. 



2.4 Numerical considerations 



While ( 42 l presents a method for numerical solution of the transition probabilities 



Pm,n(t) for a BDP with arbitrary birth and death rates, it is not yet an algorithm 
for reliable evaluation of these probabilities. In order to develop a reliable numerical 
method, we must: 1) characterize the error introduced by discretization of the integral 



in (41 1; 2) determine a suitable method to evaluate this nearly alternating sum while 
controlling the error; and 3) accurately and rapidly evaluate the infinite continued 



fraction in ( 15 1 



Abate and Whitt show that the discretization error that arises in ( 42 I is 

oo 

e d = Yl e ~ kAp ™,n {(2k + l)t) , 



(43) 



fc=l 



and when P m ,n(t) < 1, 



-kA 



,-A 



(44) 



k=l 



when e is small. Then to obtain < 10 7 , we set A = 7log(10). As 



Abate and 



Whitt point out, the terms of the series (42 1 alternate in sign when 



Re /, 



A + 2kni 
2t 



(45) 



15 



has constant sign. This suggests that a series acceleration method may be helpful in 
keeping the terms of the sum manageable and avoiding roundoff error due to summands 



of alternating sign. We opt to use the Levin transform for this purpose (Levin 1973 
Press|2007| [Numerical Recipes Software|2007 l. 



Evaluation of rational approximations to continued fractions by repeated applica- 
tion of Lemma [I] is appealing, but suffers from roundoff error when denominators are 
small (iPress 2007). To evaluate the infinite continued fraction in the summand of (42 1, 



we use the modified Lentz method (Lentz 1976 Thompson and Barnett 1986 Press 



20071. To demonstrate, suppose we wish to approximate the value of /o,o( s ), given by 



Q by truncating at depth k. Then 



0») 



B k (s) 



(46) 



is the fcth rational approximant to the infinite continued fraction /o,o( s )- I n t ne modi- 
fied Lentz method, we stabilize the computation by finding the ratios 



Cfc = ~ A and D k = — — 

A k-l £>k 



(k) 

so that /q can be found iteratively by 



(47) 



/o,0 — fo,0 ^CkDk- 



(48) 



Using Lemma [T] we can iteratively compute Cj. and via the updates 



C k = h + 



dk 



and Dk 



b k + a k D k -\ 



(49) 



In practice, we must evaluate the continued fraction to only a finite depth, but 
we must evaluate to a depth sufficient to control the error. Suppose we wish to eval- 
uate the infinite continued fraction /o,o(s) given by Q at some complex number s. 
Intuitively, we wish to terminate the Lentz algorithm when the difference between suc- 
cessive convergents is small. However, it is not immediately clear how the difference 



16 



between convergents (s) — f^L (s) is related to the absolute error ,fo.o( s ) — /no • 



Craviotto et al ( 1993 1 make this relationship clear by furnishing an a posteriori trun- 
cation error bound for Jacobi fractions of the same form as Q in this paper. Assuming 

(k) 

that /g o (s) = A k (s)/B k (s) converges to /o,o(s) as k -> 
the bound 



Craviotto et al 



(19931 give 



kow-/8>w| < 





Bk(a) 






Im 


' B fc ( s ) \| 





l/oSw-ZoSr^wl. ( 5 °) 



that is valid when Im(s) is nonzero. Note that B k {s) / B k _i{s) = 1/D k (s), so (501 
is easy to evaluate during iteration under the Lentz algorithm. Therefore, we stop at 
depth k in the Lentz algorithm when 



|Im(l/D fe ( S ))| 



(51) 



is small. 



2.5 Numerical results 



Although our error-controlled method is designed to be used when an analytic solution 
cannot be found, we seek to validate our numerical results by comparison to available 
analytic and numerical solutions. For the simple BDP with An = nX and \i n = np,, our 
numerical results agree with the values from the well-known closed-form solution given 



explicitly in Bailey ( 1964 1 as 



min(m,n) 
Prn,n(t) = ^ ^ 



3=0 



P m ,o(t) = a T ' 



m — 1 



rn\jm + n j 1 1 , a-j ^n-j q _ q _ 



(52) 



where 



(i (>-")' -l) X U*-^ - 1 

— w ~ and P = 7t r; 



(53) 
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n 

Fig. 1 Comparison of transition probabilities Pip, n ( t = 1) computed by our error-controlled 
method and that of |Murphy and O'Donohoe] | |1975[ l for the immigration-death model with 
A„ = 0.2 and fi n = 0.4n. 'i'he open circles ar e the values given by our method. The solid line 
corresponds with the approximant method of |Murphy and Q'Donohoe| with k = 2 (solid line), 
k = 3 (dashed line), and k = 4 (dotted line). In our experience, the approximant method fails 
whenever n + m + k is greater than approximately 20. It is interesting to note that increasing 
the depth of truncation k in the approximant method actually worsens the approximation. 



Murphy and O'Donohoe ( 1975 1 give numerical probabilities for four general birth- 
death models: a) immigration-death with A rl = 0.2 and fi n = 0.4n; b) immigration- 
emigration with An = 0.3, fio = 0, and fin = 0.1; c) queue with \ n — 0.6, /in = 0, 
Hi = H2 = 0.2, /i3 = Hi = 0.4, and [in = 0.6 for n > 5; and d) A n = 0.4, fin = 0.1y/n. 
Our results agree with those computed by Murphy and O'Donohoe for each of the four 
models given in Tables 2 through 7 in their paper ( |Murphy and Q'Donohoe|1975[ ). We 
note that Murphy and O'Donohoe] did not report probabilities for m > 2 or n > 5 
in any of their four models. In our experience, their method performs poorly when 
n + 7Ti + k is greater than approximately 20. 

As a demonstration of the instability of the approximant method, we contrast the 
numerical results given by our error-controlled method with those obtained using the 
approximant method, that we implemented as described in |Murphy and O'Donohoe] 



(19751, except for some rescaling of intermediate quantities to avoid obvious sources 



of roundoff error. Figure [T] shows this comparison, using model (a) above, for three 
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values of the truncation index k. Note that increasing the truncation depth k in the 
approximant method does not improve the error. 

3 Applications 

Drawing on the robustness and generality of our error-controlled method, we conclude 
with four models in ecology, genetics, and evolution whose analytic solutions remain 
elusive and where past numerical approaches have fallen short. Using our approach, 
computation of transition probabilities is straightforward, and the techniques outlined 
above may be used without modification. Some of the examples are well-known models, 
and others are novel. In some cases, the orthogonal polynomials satisfying |3| are 
known, and hence a solution could be numerically computed using pj, provided there 
are good ways of evaluating the polynomials. Often, a severe drawback of using known 
orthogonal polynomials to compute a solution based on (J2J is that the polynomials 
are model-specific. This makes experimentation and model selection difficult, since 
computation of transition probabilities depends on a priori analytic information about 
the polynomials and measure associated with the BDP. Our method does not rely on 
a priori information about the process, other than the birth and death rates for each 
state. 

3.1 Immigration and emigration 

Consider a population model for the number of organisms in an area, and suppose new 
immigrants arrive at rate v, and emigrants leave at rate 7. Organisms living in the area 
reproduce with per-capita birth rate A and die with rate \l. Define the linear rates 

An = n A + v and fi n — n/i + 7. (54) 



For the case 7 = 0, an analytic expression for the orthogonal polynomials is known 
(Karlin and McGregor 1958a I. For nonzero 7, orthogonal polynomials are available 
from which a solution of the form B may be computed (Karlin and McGregor|1958a 
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Fig. 2 Transition probabilities for the immigration/emigration model with A = 0.5, v = 0.2, 
fj, = 0.3, and 7 = 0.1. The top panel shows Pio,n(i) with t = 1 (solid line), t = 2 (dashed 
line), t = 3 (dotted line), and t = 4 (dash-dotted line) for n = 0, . . . , 50. The bottom panel 
shows Pio,n(t) with n = 15 (solid line), n = 20 (dashed line), n = 25 (dotted line), and n = 30 
(dash-dotted line) for t e (0,20). 



Ismail et al|1988 l. However, using our error-controlled method, we can easily find the 
transition probabilities without additional analytic information. Figure [2] shows an 
example of the time-evolution of Pio,n{t) for various times t and states n, with the 
parameters A = 0.5, v — 0.2, /i = 0.3, and 7 = 0.1. The approximant method method 



of Murphy and O'Donohoe fails to produce useful probabilities for n > 10 (not shown). 
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Allee Effect Logistic growth 



Logistic decay 




Fig. 3 Behavior of logistic/ Allee model. The upper panel shows a plot of birth (solid line) 
and death (dashed line) rates for states n = 0, . . . , 60, and parameters \ = 1, /j, = 0.1, M = 20, 
a = 0.2, and fj = 0.3. The different phases of growth are labeled in the shaded regions. The 
lower panel shows stochastic realizations of the logistic/ Allee model for various starting values. 
The shaded regions correspond with the shaded phases of growth in the upper panel. 



3.2 Logistic growth with Allee effects 



Populations of organisms that occupy a finite space may be subject to various con- 
straints on their growth. The per-capita birth rate may decline when there are more 
organisms than the ecosystem can sustain (Tan and Piantadosi 19911. This can hap- 
pen when there are too many organisms competing for the same food supply. The 
decay of population size above some carrying capacity is usually called logistic growth 
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Fig. 4 Logistic/ Allee model probabilities of extinction P m ,a(t) for initial population sizes 
m = 1 (solid line), m = 5 (dashed line), m = 10 (dotted line), and m = 15 (dash-dotted line). 
The full model parametrization is found in the text. 



by ecologists. Another density-dependent constraint is known as the Allee effect, in 
which per-capita birth rate increases superlinearly with n once a small population has 
been established, due to favorable consequences of density, such as cooperation and 
mutual protection from predators ( Allee et al| [T949). As a realistic example of a gen- 
eral BDP that has no obvious solution by orthogonal polynomials, we seek a model 
that both transiently supports growth above the carrying capacity, and reflects these 
two density-dependent constraints, similar in spirit to models described by |Tan and] 



Piantadosi (1991) and Dennis (2002) 



Qualitatively, if the per-capita birth rate with no density effects is A, then the total 
birth rate should rise faster than nX when n is small, slower than nX for intermediate n 
near the carrying capacity, and should decay toward zero for n greater than the carrying 



capacity. Tan and Piantadosi introduce a logistic birth rate X n = nX (l 



N 



for a finite 



state space model that takes values {0,1,..., N}. However, to allow for temporary 
growth beyond the carrying capacity, we choose X n oc Xn 2 e~ an for intermediate and 
large n. To achieve attenuated growth for small n as well, we scale this rate by a logistic 
function, yielding 



A n 



Xn z e' an 

I + e l3(n-M) 



and [in = nfi, 



(55) 
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where M is the population size with highest birth rate, and the death rate is assumed to 
be proportional only to the number of existing individuals. Figure[3]shows the resulting 
rates for various states n, with the different phases of population change shaded. To 
illustrate that the model produces the desired behavior, several realizations of the 
process are given in the lower panel for various starting values. The shaded regions 
correspond with the three phases of growth. Note that most paths in the lower panel of 
Figure [3] center near n — 27, where the birth rate and death rate are equal. The lower 
panel corresponds with Figure 1 in Dennis (2002). Figure [4] demonstrates the success 
of the error-controlled method in computing time-dependent extinction probabilities 
Pm,o f° r various starting values with A = 1, a = 0.2, /3 = 0.3, M = 20, and /j, = 0.1. 



3.3 Moran models with mutation and selection 

The probability of fixation or extinction of an allele in finite populations is frequently of 
interest to researchers in genetics. However, publications often rely on the probability 
of eventual extinction P m ${t — > 00), or the probability of fixation of a novel mutation 
in a population of constant size N, Pi pf(t — > 00). While these asymptotic probabil- 
ities do reveal important properties of the underlying models, the information they 
provide about the distribution of time to fixation/extinction is incomplete. In practice, 
researchers may observe that m organisms in a sample exhibit a certain trait at a 
certain time. Then P m $(t), the probability of extinction of that trait at finite times t 
in the future should presumably be of great interest, since researchers cannot reliably 
observe the process for infinitely long times. Additionally, the finite-time probability 
of fixation/extinction may exhibit threshold effects or unexpected dynamics that are 
not revealed by the asymptotic probability of such an event. 



Moran ( 1958 1 introduces a model for the time-evolution of a biallelic locus when the 
population size is constant through time. A biallelic locus is a location in an organism's 
genome in which two different genetic variants or alleles exist in a population. We are 
interested in how the number of individuals carrying each allele changes from generation 



to generation. Krone and Neuhauser ( 19971 exploit the Moran model to derive a BDP 
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counting the number of individuals with a certain allele in the context of ancestral 
genealogy reconstruction in which one allele offers a selective advantage to individuals 
that carry it. Selection greatly complicates the problem and remains an active area of 
research. In a limiting case, this process corresponds to Kingman's coalescent process 
when there is no mutation or selection (Kingman I98 2a|b I. 



To construct the Moran process with mutation and selection, suppose a finite pop- 
ulation of N haploid organisms has 2 alleles at a certain locus: Ai and A2 . Individuals 
that carry Ai reproduce at rate a and A2 individuals reproduce at rate j3. Suppose 
further that individuals carrying the A\ allele have a selective advantage over individ- 
uals carrying A2, so a > /3. When an individual dies, it is replaced by the offspring of a 
random parent chosen from all N individuals, including the one that dies. This parent 
contributes a gamete carrying its allele that is also subject to mutation. Mutation from 
A\ to A2 happens with probability u and in reverse with rate v. The new offspring 
receives the possibly mutated haplotype and the process continues. 

Let X{t) be a BDP counting the number of A\ individuals on the state space 
n £ {0, . . . ,N}. To construct the transition rates of the process, suppose there are 
currently n individuals of type A\. We first consider the addition of a new individual 
of type A\, so that n — > n + 1. For this to happen, the individual that dies must be 
of type A2- If the parent of the replacement is one of the n of type Ai, the parent 
contributes its allele without mutation, and this happens with probability I — it. If the 
parent of the replacement is one of the N — n of type A2, the parent contributes its 
allele, which then mutates with probability v. Therefore, the total rate of addition is 



_ N-n 



n . _ N — n 



(56) 



for n — 0, . . . , N with \ n = when n > N. Likewise, the removal of an individual of 
type A\ can happen when one of the n individuals of type A\ is chosen for replacement. 
If the parent of the replacement is one of the iV — n of type A2 , the parent contributes 
A2 without mutation, with probability 1 — v. If the parent is one of the n of type A\, 
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the allele must mutate to A2 with probability u. The total rate of removals becomes 



Mn = 



. n 



(57) 



for n = 1, . . . , N with = Ajy = and fi n — when n > N. Note that if n > 0, then 
Ao > so the allele cannot go extinct. Also, if u > 0, then /ijy > 0, so the Ai allele 
cannot be fixed in the population. 



Karlin and McGregor ( 1962 1 derive the relevant polynomials and measure for the 



Moran process described above, but without selection, so that a — f3. Donnelly ( 1984 \ 



gives expressions for the transition probabilities in the case where a = = 1, noting 
that when selection is introduced (via differing a and /3), his approach is no longer 
fruitful. Using our technique, computation of the transition probabilities under selection 
is straightforward. The upper panel of Figure [5] shows the probability of fixation by 
time t. The lower panel shows the finite-time fixation probability of A\, P m ,lOo{t), with 
u — so the state n = 100 is absorbing. 

Since the state space in the Moran model is finite, it is natural to consider the 
matrix exponentiation method discussed in the Introduction. We write the stochastic 
transition matrix as 



-Ao Ao 

Mi -(Ai + Mi) Ai 

M2 -(A2 + M2) A 2 



\ 



(58) 



-(Ajv + MaO Aat , 



where X n and fi n are defined by (56 1 and 1 57 1, respectively. In our experience, the 



matrix exponentiation method often works well, and its computational cost is similar 
to that of our error-controlled method. However, it is highly sensitive to rate matrix 
conditioning. For example, Figure [6] shows a comparison of transition probabilities 
from the error-controlled method and the matrix exponentiation method for the Moran 
model with N = 100, a = 210, j3 = 20, u — 0.002, and v = 0. In evolutionary terms, 
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Fig. 5 Transition probabilities for the Moran model with selection. The upper panel shows 
the probability of n individuals having allele A\ at time t, Pso,n(t) for the Moran model with 
N = 100, starting from m = 50 with u = 0.02, v = 0.01, a = 60, and = 10. We show the 
probabilities for t = 1 (solid line), t = 3 (dashed line), t = 5 (dotted line), t = 8 (dash-dotted 
line). Note that although the states and 100 are not absorbing, the mutation rates u and 
v are small enough that probability accumulates significantly in these end states. Note also 
the asymmetry in the distribution at longer times. The lower panel reports the probability 
of fixation by time t, Pm.iooWi f° r the same model, but with u = so the state n = 100 is 
absorbing. The probabilities shown are for m = 70 (solid line), m = 50 (dashed line), m = 20 
(dotted line), and m = 1 (dash-dotted line). Note the starkly different time-dynamics for 
different starting values. 



this means that mutation from A\ to A2 is impossible, and the A2 haplotype suffers 
from low fitness. Computationally, this has the effect of making fi n small for most n, 
and hence the rate matrix grows ill-conditioned. 
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Fig. 6 Comparison of Moran model transition probabilities Pso,n(i = 0.2) computed by 
two methods with N = 100, a = 210, j3 = 20, u = 0.002, and' v = 0. The open circles 
correspond with our error-controlled method, and the solid line corresponds with the matrix 
exponentiation method. This choice of parameters causes wild fluctuations in probabilities 
reported by the matrix exponentiation method since the stochastic rate matrix becomes nearly 
singular. 



Although the rate matrix in this example is nearly defective, this choice of parame- 
ter values is not unreasonably extreme. For example, researchers in population genetics 
often wish to test the hypothesis that selection occurs in a dataset. They fit parame- 
ters for models with selection (full model) and without selection (restricted model) and 
perform a likelihood ratio test of this hypothesis. If the estimates of /3 and u in the full 
model are small, they may be unable to reliably compute the probability (likelihood) 
of the data, given the estimated parameter values under the full model. 



3.4 A frameshift-aware indel model 



Thorne et al ( 1991 1 introduce a BDP modeling insertion and deletion of nucleotides 
in DNA for applications in molecular evolution. The authors model the process of 
sequence length evolution by assuming that a new nucleotide can be inserted adjacent 
to every existing nucleotide, and every existing nucleotide is subject to deletion, at a 
constant per-nucleotide rate. This corresponds to the simple BDP with \ n — n\ and 
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fin = nfi. If a sequence has m nucleotides at time and there are n nucleotides at time 
t later, the probability of this event is P m ,n(t)- 

However, an important aspect of biological sequence evolution is conservation of 
the structure and biophysical properties of proteins that result from transcription and 
translation of DNA sequences. After coding DNA is transcribed into RNA, ribosomes 
translate 3-nucleotide chunks (codons) of the RNA into a single amino acid residue, 
that is then joined to the end of a growing protein polymer. Insertions or deletions 
(indels) in a DNA sequence that result in a shift in this triplet code are called "frame- 
shift" mutations. It is likely that a frame-shift indel occurring in a protein-coding DNA 
sequence results in a protein that is prematurely terminated or possesses structural 
and chemical characteristics unlike the ancestral protein. Insertions or deletions whose 
length is a multiple of three should be more common. We seek to model this behavior 
in a novel way: suppose the indel process is a BDP similar in spirit to the one presented 



by Thorne et al (19911, and the rate of insertion and deletion of nucleotides depends 



on the number of nucleotides already inserted, modulo (mod) 3: 



A„ = < 



n/?o if n — 1 = mod 3 
nfii if n - 1 = 1 mod 3 and M« = ' 
n/?2 if n — 1 = 2 mod 3 



njQ if n — 1 = mod 3 

nji if n — 1 = 1 mod 3 • (59) 

nj2 if n — 1 = 2 mod 3 



Here we assume that fa > /3o,f3%, and 71 > 70,72 so that transitions to state n 
such that n — 1 = mod 3 occur at a faster rate per nucleotide. The linear-periodic 
nature of these birth and death rates make solution of the orthogonal polynomials and 
measure corresponding with this BDP difficult. The approximant method of |Murphy] 
|and Q'Donohoe| also fails here for large n. However, using our error-controlled method, 
numerical results are readily available. Figure [7] shows P\ t n{t) for n = 0, . . . , 50 at 
various times t. Note that the distribution of the number of inserted bases has peaks 
at the integers mod three. Finally, it is worth noting that the dearth of tractable BDPs 
for indel events has been a major deterrent in statistical sequence alignment and we 
are actively exploring solutions to this problem using our error-controlled method. 
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Fig. 7 Framcshift-aware indcl model probability of observing n inserted DNA bases, given 
starting at m = 1. The transition probability Pi, n (t) is shown for t = 5 (solid line), t = 7 
(dashed line), t = 9 (dotted line), and t = 11 (dash-dotted line), with parameters ft) = 0.3, 
ft = 1, ft = 4, 70 = 2, 7i = 0.2, and 72 = 0.2. 



4 Conclusion 



Traditionally the simple BDP with linear rates has dominated modeling applications, 
since its transition probabilities and other quantities of interest find analytic expres- 
sions. However, increasingly sophisticated models in ecology, genetics, and evolution, 
among other fields, may necessitate more advanced computational methods to handle 
processes whose birth and death rates do not easily yield analytic solutions. We have 
demonstrated a flexible method for finding transition probabilities of general BDPs 
that works for arbitrary sets of birth and death rates {A n } and {/in}, and does not 
require additional analytic information. This should prove useful for rapid develop- 
ment and testing of new models in applications. For simple models whose solution is 
available, we find that our method agrees with known solutions and remains robust 
for large starting and ending states and long times t. It is our hope that the method 
presented here will assist researchers in understanding the properties of increasingly 
rich and realistic models. 
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5 Appendix 



5.1 Approximant method 



Murphy and O'Donohoe ( 1975 1 approximate the inverse Laplace transform of ( 15 1 
by first truncating the continued fraction as a rational approximant through a partial 
fractions sum. To illustrate the pitfalls of this approach, we derive the inversion expres- 
sions presented by |Murphy and Q'Donohoe| and analyze their properties. We provide 
an example to show that this technique can become numerically unstable. We first seek 
to uncover the truncation error in the time domain of the transition probabilities. If 



we truncate the continued fractions (15 1 at depth k, we have 



fL k Us)=l J] N \ *J nam + 2 ^ b -r^ forn<m,and 

I - LJ - ' B m+ i+ b m+ 2+ m +3+ rn+k 



n-1 



(60) 



A k ) t \ I TT \ I B m B n a n+ 2 dn+3 a n+k c ^ 

fm n(s) = II — • ■ • - for n > m 

' ' B n+1 + b n+2 + b n+3 + 6„.i. 



+k 



For concreteness, suppose in what follows that n > m. Note that the denominator 

(n) 

of the second equation is simply B n+k . Let Al be the numerator of the continued 



fraction in the second equation in (601, so 




f^n = I II ^ I (61) 



where A k satisfies A k = A k , A ™ — Y\!j=\ a ji an( i 



A ( k n) = a n+k A k % + b n+k A k % (62) 
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Note also that the difference between truncated estimates in the Laplace domain (s) is 

An+k A n _ A n+k B n ~ A n B n+k 



Bn+k B n B n+k B n 

= (-l)"4"> _ 

B n +kB n 



This yields the generalized determinant formula 



(63) 



A n+k B n - A n B n+k = (64) 



and at a root Sj of B n+k (s), we have 

A k n \ Si ) = (-l) n A n+k ( Sl )B n (sA. (65) 
Now if si, S2, ■ ■ ■ , s n are the roots of B n (s), we have, using the previous line and a 



partial fractions decomposition of (601, the formula for the Laplace transform of the 



transition probability P m ,n(t), truncated at k, 



\J=m 



B n +k(s) 



11 Aj ' jgM (66) 



' n— 1 \ n+/e 



rr A I y -B m (s)-Bn(si)^n + fc (si) 



since we only require the values of A n+k (s) and B n (s) at the zeros of B n+k (s). Then 
inverse transforming, an approximate formula for the transition probability P m ,n(t) is 

m ~ I TT X I V" B m(Si)Bn{Si)A n+k {Si) s . t 

m ' ni) ~\iL s h • (6T) 
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The roots of B n (s), used in (661 and |67|, are often found numerically as follows. 
Consider the characteristic polynomial det (_B n + s/) of the matrix 



Br, 



Ao 1 
Aomi Ai 1 

Ai/x 2 A 2 + ^2 



V 



A„-3A*n-2 A„_2 + 1 



(68) 



It is clear that the nth partial denominator B n (s) — det(B n + si), and this quantity 
is zero when —s is an eigenvalue of the matrix B ri . Therefore, the negatives of the 
eigenvalues of B n are the roots of B n (s). Furthermore, B n can be transformed into 
a real symmetric matrix via a similarity transform and hence B n (s) has precisely n 
roots, all of which are simple, real, and negative. One usually finds these eigenvalues via 
the QR algorithm or similar numerical techniques (Press 2007). However, the iterative 



eigendecomposition of (68 1 generates small errors in the eigenvalues for large n + k. 



These errors are amplified in the product in the denominator of each summand in (671 



resulting in a sum with both positive and negative terms that may be very large. !Klar 



et al (20101 encounter similar instability in this algorithm. Their solution is to find 



the roots of the terms in the numerator and compute each summand as a product 
of individual numerators and denominators in an attempt to keep roundoff error in 
the product from accumulating. So, if z\, . . . , z n+ f, are the roots of A n+ ^ then ( |67[ ) 
becomes 



(n— 1 \ n+fe 
II X l B m{si)B n { Si ){ Zi - Si ) Yl 

j=m I i=l jjt% 



-Sit 



(69) 



This procedure does improve the numerical stability of the computation, but requires 
two eigendecompositions of possibly large matrices for every evaluation of P m ,n(t), 
increasing the computational cost and, for large m and n, the roundoff error. In our 



:',2 



opinion, it is more advantageous to avoid truncation of the continued fraction |9| at 
a pre-specified index, and instead evaluate the continued fraction until convergence 
during numerical inversion. Figure [I] shows how approximant methods fail for large n. 



5.2 A power series method 



Parthasarathy and Sudhcsh (2006a) present exact solutions by transforming continued 



fractions such as |9| into an equivalent power series. Wall ( 1948 1 shows that Jacobi 
fractions of this type can always be represented by an equivalent power series. However, 
the small radius of convergence of power series expressions for transition probabilities 
can limit their usefulness for long times or large birth or death rates. Parthasarathy 



and Sudhesh show that Po, n (t) has a power series representation given by 

a—l \ oo 
Pm,n{t) = | II ) E(-l) m %,2™) ( 



k=0 



m=0 



where 



il + l 



i 2 + l 



A(m,n) = ^ a n ^ a. h ^ 



(m + n)\ ' 



im-l + l 



(70) 



(71) 



i 1= i 2 =0 i 3 =0 i m =0 

with A(0, n) = 1 ( jParthasarathy and S udhcsh 2006a bj. Here, ai n — \ n and 02n+l = 
jj, n in the notation used in their papers. This approach is unique because it yields an 
exact analytic expression for the transition probabilities of a general BDP. However, the 
radius of convergence of the power series depends on the specified rates, and this radius 
may be quite small. To illustrate the pitfalls of this approach, consider a n = (n + 1)A, 
corresponding to the BDP with An = (2n + 1)A and /i„ = 2nA ( jParthasarathy and| 



Sudhesh 2006a Example 4.6). The power series for the transition probability in this 
process becomes 



P ,n(t) = E 



n-\-m 



. (2n + 2m)! (Xt/2) 
m\n\ (n + m) 



(72) 
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Then the radius of convergence R of the power series is given by 

. n+m+l 



1/R 



lim 



(2n + 2m + 2)! (§) 



X 



m\n\(n + m) 



— lim 

m— >oo 



oo (m + l)!n!(n + m+l)! (2n + 2m)!(§) 

(2m + 2n + l)(2n + 2m + 2) / A 
(m + l)(n + m + l) \2 



n+m 



2m + 2n + 1 , 
lim A 

m— s-oo 771 + 1 

2A. 



(73) 



And so the series diverges when 2Xt > 1. To illustrate the limitations of the power 
series approach, note that in this process, the transition intensity from to 1 is A, so 
the expected first-passage time from to 1 is E(Xn.i) = 1/A. Therefore, we cannot 



evaluate (72l when t is greater than E(Tq i)/2. If n is much greater than 1, we may be 



unable to reliably evaluate Po,n(£) for times near E(Tn in ) 
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